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ABSTRACT:  A  computer  program  Is  presented  for  calculating  the 
various  components  comprising  the  total  drag  coefficient  for  slightly 
blunted  slender  cones.  The  program  is  applicable  for  either  the  case 
of  a  completely  laminar  boundary  layer  or  where  transition  occurs 
on  the  cone.  Weak  viscous  interaction  effects  are  taken  into  account 
Local  flow  properties  which  exist  at  the  outer  edge  of  the  boundary 
layer  along  the  cone  are  also  computed. 
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INTRODUCTION 


A  program  has  been  formulated  and  coded  for  use  on  an  1 W  7090 
for  calculating  the  drag  coefficient  of  a  slightly  blunted  cone. 

The  program  is  applicable  for  cases  where  the  boundary  layer  is 
either  completely  laminar,  completely  turbulent,  or  the  flow  is  mixed, 
with  transition  occurring  on  the  cone. 

The  equations  that  are  used  do  not  require  that  the  flow  along 
the  outer  edge  of  the  boundary  layer  be  isen tropic.  A  variation  in 
entropy  along  this  surface  due  to  the  curved  bow  shock  wave  is 
allowed.  This  effect  results  in  a  decrease  in  the  skin-friction 
drag  from  that  for  a  sharp  cone.  For  slender  cones  where  the  skin- 
friction  drag  is  a  significant  portion  of  the  total  drag,  this  effect 
of  blunting  the  cone  plays  an  important  role  in  the  determination  of 
the  total  drag. 

Viscous  interaction  effects  are  considered,  and  the  increase  in 
skin-friction  coefficients  due  to  transverse  curvature  and  induced 
pressure  effects  are  taken  into  account.  The  increase  in  pressure 
drag  due  to  Induced  pressure  resulting  from  large  boundary-layer 
displacement  thicknesses  is  also  taken  into  account.  These  effects 
can  be  significant  for  conditions  of  low  Reynolds  numbers  and  high 
Mach  numbers. 

As  a  result  of  the  method  used  to  calculate  the  skin-friction 
coefficients,  local  flow  properties  at  the  outer  edge  of  the  boundary 
layer  are  also  determined.  The  integrated  boundary-layer  thickness 
parameters  (momentum  and  displacement)  are  calculated. 

No  attempt  was  made  to  Include  real  gas  effects  such  as  dis¬ 
sociation  or  ionization  of  the  flow  over  the  cone. 


The  program  can  be  used  to  evaluate  the  effects  of  Mach  number, 
Reynolds  number,  wall-to-recovery  temperature  ratio,  nose  bluntness, 
cone  angle,  and  the  location  of  boundary-layer  transition  on  the 
total  drag  coefficient  or  on  the  various  components.  Some  of  these 
effects  will  be  shewn  in  comparisons  of  calculations  with  experi¬ 
mental  data,  and  in  discussions  of  the  drag  characteristics  of  slightly 
blunted  slender  cones. 


Care  has  been  taken  in  coding  the  program  for  the  machine  cal¬ 
culations  to  provide  as  simple  an  input  as  possible  to  afford  ease 
in  operation.  Information  concerning  running  of  the  program  is  pro¬ 
vided  in  Appendix  A,  along  with  a  complete  listing  of  the  program. 
Appendix  B  is  a  print-out  of  a  sample  test  case  which  was  run  using 
this  program. 


DESCRIPTION  OF  EQUATIONS 

The  total  drag  coefficient  is  formed  by  summing  the  pressure 
drag  coefficient,  the  skin-friction  drag  coefficient,  and  the  base 
drag  coefficient  plus  Incremental  increases  to  the  pressure  and 
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akin-friction  drag  coefficients  due  to  viscous  interaction  effects. 

Hie  equations  used  to  obtain  each  of  these  coefficients  will  be 
discussed  separately.  As  will  be  seen,  to  perform  some  of  these 
calculations  it  is  necessary  to  have  values  of  inviscid  sharp-cone 
properties  for  the  particular  case  of  interest.  The  equations  used 
for  these  computations  will  also  be  discussed. 

Inviscid  Cone  Properties 

Hie  following  method  is  used  cu  calculate  the  flow  properties 
which  would  exist  on  a  hypothetical  sharp  cone  having  the  same  half 
angle  as  the  blunted  cone  being  considered,  and  for  the  same  free- 
stream  Mach  number  being  considered.  These  properties  are  calculated 
using  the  simple  relations  given  by  Blick  in  reference  (3).  Relations 
are  presented  for  calculating  the  temperature  and  the  velocity  of  the 
inviscid  flow  along  the  cone  referenced  to  the  respective  free-stream 
quantities.  Using  these  two  relations,  the  Mach  number  can  be  deter¬ 
mined.  The  inviscid  sharp-cone  pressure  can  be  calculated  using  the 
relation  given  by  Llnnell  and  Bailey  (ref.  (4))  and  presented  again 
in  reference  (3).  These  equations  are  used  in  the  program  for  cal¬ 
culating  the  inviscid  cone  pressure,  temperature,  and  Mach  number. 

Geometric  Characteristics  of  a  Spherically  Blunted  Cone 

It  was  desired  to  be  able  to  express  geometric  quantities  for 
the  spherically  blunted  cone  in  terms  of  the  cone  half  angle,  the 
nose-to-base  radius  ratio,  and  the  radius  of  the  base.  Two  quanti¬ 
ties  in  particular  are  necessary:  the  distance  from  the  stagnation 
point  to  the  point  of  tangency  between  the  spherical  nose  cap  and  the 
conical  frustum,  Sqi,  and  the  total  wetted  length  from  the  stagnation 
point  to  the  base  of  the  cone,  Sj,^.  These  two  quantities  are  given  by 


sT  -  rB  fN  /n  .  e  ) 
rB  V8  j 

and 

SMAX  "  rB 


rN  tt  _  0 „  j  1  **  rN/rB  SIN(i  ~  8c) 

55  2  - sinV - 


Laminar  Boundary-Layer  Calculations 

Calculations  of  the  laminar  boundary-layer  characteristics  are 
first  made  neglecting  viscous  interaction  effects  such  as  induced 
pressure  and  transverse  curvature.  The  equations  used  for  making 
these  laminar  calculations  were  derived  by  Wilson  and  reported  in 
reference  (5).  This  is  a  momentum- integral  method  which  utilized 
available  flat  plate  results.  The  flow  characteristics  are  followed 
along  streamlines  through  the  bow  shock  wave,  using  oblique  shock 
relations,  and  then  expanded  lsen tropically  to  the  point  where  the 
streamline  intersects  the  outer  edge  of  the  boundary  layer.  The 
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boundary-layer  thickness  defining  the  outer  edge  is  determined  through 
a  mass  balance  of  the  flow  in  the  free  stream  and  that  in  the  boundary 
layer  at  this  intersection  point.  Having  defined  the  properties  at 
the  outer  edge  of  the  boundary  layer,  the  boundary-layer  character¬ 
istics  are  determined  from  a  step-by-step  technique  using  the  momentum- 
integral  equations  and  necessary  thickness  parameters  and  skin-friction 
coefficients  taken  from  flat  plate  results.  The  viscosity  ratios  were 
evaluated  using  Southerland  * s  equation  in  the  form 

u  .  1  +  198.6/%  fT  \3/a 

hi  T/T].  +  19«.6/T1 


The  momentum- thickness  Reynolds  number  based  on  flow  properties  Just 
outside  the  boundary  layer  is  calculated  using  the  relation 


Re 


1 


Rec  «1  (%  +  198.6)  /tc\2 


These  laminar  equations  are  applied  to  the  region  between  the 
point  of  intersection  of  the  spherical  nose  and  the  conical  frustum 
and  either  the  end  of  the  cone,  in  the  case  of  an  all-laminar  boundary 
layer,  or  the  point  of  boundary-layer  transition  for  the  case  of  a 
mixed  laminar- turbulent  boundary  layer. 


Turbulent  Boundary-Layer  Calculations 

For  cases  where  transition  from  laminar  to  turbulent  flow  is 
specified  to  occur  on  the  body,  then  -at  this  point  of  transition  the 
calculation  procedure  is  changed  to  one  applicable  to  the  turbulent 
boundary  layer.  Sie  procedure  used  is  still  a  momentum-integral 
method,  and  was  formulated  and  reported  by  Wilson  in  reference  (6). 
Having  values  for  Re  ,  RXq,  and  Cf  computed  at  the  point  of  transition 

using  the  laminar  boundary-layer  equations  just  described,  a  value  for 
Rq  ,  some  distance,  x,  downstream  from  the  transition  point  can  be 

calculated  using  an  integral  form  of  the  momentum-integral  equation 
as  described  in  reference  (6).  This  then  allows  a  step-by-step, 
simultaneous  solution  of  the  turbulent  boundary-layer  characteristics 
and  the  lnviscid  flow  field  properties  at  the  outer  edge  of  the 
boundary  layer  to  be  carried  out.  As  indicated  in  reference  (6),  to 
obtain  the  thickness  parameters  it  is  necessary  to  assume  a  velocity 
profile.  A  power  law  profile  is  assumed  with  an  exponent  of  1/7, 
although  it  is  a  simple  matter  to  perform  the  calculations  using  a 
different  value  for  the  exponent.  Some  errors  were  found  in  one  of 
the  equations  in  reference  (6).  This  has  been  corrected  in  this 
program,  and  the  correct  form  of  the  equation  appears  in  the  listing 
of  equations  in  a  following  section  of  this  report. 
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Viscous  Interaction  Effects 


At  high  Mach  numbers  and  low  Reynolds  numbers  the  thickness  of 
a  laminar  boundary  layer  on  a  slender  cone  can  be  significant  com¬ 
pared  to  the  local  radius  of  the  cone.  Under  these  conditions  an 
interaction  between  the  thick  boundary  layer  and  the  Inviscid  flow 
field  can  exist.  Two  predominate  viscous  interaction  effects  are 
induced  pressure  and  transverse  curvature.  Ihese  viscous  interaction 
effects  are  characterized  by  a  viscous  interaction  parameter, 
defined  as 


where 


X 


M5  -fc" 

“f'Te 


C  « 


T  +  198.6  -jTw/T 

'V+T98T6 - 


For  small  values  of  Y  {weak  interaction),  Probstein  (ref.  (4))  has 
presented  equations  for  calculating  the  induced  pressure  on  a  cone. 

Also  in  reference  (4),  Probstein  has  presented  equations  for  calculating 
the  increase  in  the  local  skin-friction  coefficient  on  a  cone  due  to 
induced  pressure  and  transverse  curvature  effects.  The  relations  in 
reference  (4)  were_developed  for  a  sharp  cone,  assuming  that  the  prop¬ 
erties,  including  xc»  were  constant  along  the  cone.  Since  for  these 
present  calculations  the  Mach  number,  and  hence  x*  vary  significantly 
along  the  cone,  the  local  value  of  xe  at  points  along  the  cone  are 
used  in  the  calculations  for  the  induced  pressure  and  the  increase 
in  the  skin-friction  coefficients,  instead  of  5Cc  based  on  inviscid 
sharp-cone  values  as  shown  in  reference  (4).  The  equations  used  to 
calculate  the  induced  pressure  effects  involve  two  functions,  F-^  and 
Fg,  both  of  which  are  functions  only  of  K,  the  hypersonic  similarity 
parameter.  Bie  functions  F]_  and  F2  are  presented  in  reference  (4)  in 
both  equation  form  and  graphically.  Due  to  the  complexity,  however, 
of  the  equations,  the  curves  were  approximated  for  values  of  K  between 
1  and  2  with  a  simple  linear  expression  for  F^^  and  a  quadratic  expres¬ 
sion  for  F2.  These  two  expressions  are 

FX(K)  -  O.8992  -  0.1049KC 


and 

F2(K)  -  0.6192  -  0.3822Kc  +  0.0798KC2 

where 

kc  -  M»ec 

Drag  Coefficient  Components 

In  all  cases,  the  calculated  drag  coefficients  were  based  on  free- 
stream  dynamic  pressure  and  the  area  of  the  base  of  the  cone,  A  skin- 
friction  drag  coefficient  is  obtained  by  integrating  the  local  shear 
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stress  over  the  surface  of  the  cone  and  then  dividing  the  results  by 
the  just-mentioned  reference  quantities.  Since  the  local  skin- 
friction  coefficients  calculated  using  the  equations  given  in  refer¬ 
ences  (4),  (5),  and  (6)  are  referenced  to  local  dynamic  pressure 
and  the  exposed,  wetted  area,  then  the  shear  stress  is  simply 

t  «  Cfqe  Awetted 

The  skin-friction  drag  coefficient  can  then  be  written  as 

■S 

rcfds 

I 

where 

r  -  H  £  &  *b  +  Is  -  st)sin  ,]  cos  ec 

Bie  increase  in  the  skin-friction  drag  coefficient  due  to  induced 
pressure  is  given  by 

4CDpip  ’  k  f  r(AIV dfl 

B  sj 

The  increase  in  the  skin-friction  drag  coefficient  due  to  transverse 
curvature  is  given  by 


„  2 
°Dp  ■  r| 


The  pressure  drag  coefficient  for  the  conical  frustum  is  obtained 
by  integrating  the  longitudinal  component  of  the  local  pressure 
coefficient  over  this  portion  of  the  bcdy  and  then  dividing  the  results 
by  the  free-stream  dynamic  pressure  and  the  base  area.  Hence,  the 
pressure  drag  coefficient  can  be  written  as 

2  0s 

cDp  »  4  Cp  TSIN  9C  ds 

rB  J  e 

°T 


The  increase  in  the  pressure  drag  coefficient  due  to  induced 
pressure  is  given  by 
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V  "  ^  l  s  fr  D"  ^  +  '^J8"  e° 

For  the  spherical  nose  cap  the  pressure  is  calculated  assuming 
that  modified  Newtonian  impact  theory  is  valid.  The  pressure  coeffi¬ 
cient  is  given  as 


CP  -  °PM  ^ 


where  Cp^  is  the  maximum  pressure  coefficient  occurring  at  the 

stagnation  point  and  is  assumed  equal  to  the  stagnation  pressure 
resulting  from  bringing  the  flow  to  rest  isentroplcally  after  passing 
through  a  normal  shock.  The  angle  co  is  that  angle  between  the  free- 
stream  velocity  vector  and  a  tangent  to  the  surface  at  the  point 
being  considered.  The  resulting  pressure  is  integrated  over  this 
spherical  nose  cap  and  then  divided  by  the  free-stream  dynamic  pres¬ 
sure  and  the  base  area  of  the  cone.  This  integration  can  be  performed 
in  closed  form  resulting  in  the  relation 


A  value  for  the  base  drag  coefficient  is  calculated  assuming 
a  vacuum  exists  at  the  base,  as  usually  done  at  hypersonic  Mach 
numbers,  and  is  given  by 


COMPARISON  OF  THEORETICAL  AND  MEASURED  DRAG  COEFFICIENTS 

Calculations  have  been  performed  to  compare  theoretically 
obtained  total  drag  coefficients  with  experimentally  obtained  data. 
Most  of  the  data  were  obtained  in  the  NOL  ballistics  range  facili¬ 
ties,  although  some  wind-tunnel  data  are  also  used.  Tests  have  been 
chosen  so  that  comparisons  can  be  made  over  a  variation  of  Mach 
number,  Reynolds  number,  nose-to-base  radius  ratio,  and  wall-to- 
stagnation  temperature  ratio.  Although  there  are  data  with  some 
variation  in  cone  angle,  this  variation  is  small.  For  all  of  the 
calculations  it  was  assumed  that  the  base  pressure  was  constant 
over  the  base,  and  ea.ual  to  one  tenth  of  the  free-stream  static 
pressure.  The  agreement  between  computed  and  measured  drag  coeffi¬ 
cients  is  very  good,  with  the  exception  of  one  set  of  data  which 
will  be  discussed  later.  The  following  is  a  detailed  description 
of  all  of  the  various  comparisons. 

A  comparison  between  calculated  and  measured  drag  coefficients 
for  a  6.3-degree  half -angle  cone  is  presented  in  figure  1.  These 
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data  were  obtained  at  nominal  Mach  numbers  of  9  and  13.  Three  curves 
are  shown  in  figure  1.  The  data  at  Mach  number  9  were  presented  in 
reference  (7),  while  the  data  at  Mach  number  13  were  presented  in 
reference  (8).  The  blunted  cones  used  in  these  tests  had  half  angles 
of  6.3  degrees  and  nose-to-base  radius  ratios  of  approximately  0.03 
and  0.30.  Calculations  were  performed  for  each  set  of  data  for  the 
conditions  indicated  in  figure  1.  Data  from  tests  of  the  most  blunt 
cone  are  indicated  by  the  square  symbols.  Agreement  between  the 
calculated  values  and  the  experimental  values  is  seen  to  be  good  , 
over  a  Reynolds  number  range  from  2  x  105  to  approximately  3  x  10°. 
These  Reynolds  numbers  are  based  on  free-stream  properties  and  axial 
length  of  the  body.  Data  for  the  sharper  configuration,  tested  over 
approximately  this  same  Reynolds  number  range,  are  seen  to  be  con¬ 
sistently  higher  than  the  calculated  values  by  approximately  13  per¬ 
cent.  These  data  have  been  reviewed  with  the  author  of  reference  (7). 
It  was  determined  that  one  possible  source  of  this  deviation  could 
have  resulted  from  the  nose  bluntness  being  larger  than  reported. 

The  models  used  in  these  tests  had  a  design  nose  diameter  of  only 
0.007  inch.  This  was  extremely  difficult  to  maintain  during  manu¬ 
facturing  both  in  size  and  shape.  Further,  a  nose  of  this  size  is 
very  susceptible  to  blunting  during  flight  due  to  aerodynamic  heating. 
An  inspection  of  some  shadowgraphs  made  during  these  tests  revealed 
that  nose  diameters  could  be  as  large  as  0.020  inch.  Calculations 
have  been  made  for  a  nose-to-base  radius  ratio  of  0.1  and  are  shown  as 
the  dotted  curve  in  figure  1.  It  is  seen  that  this  curve  is  in  much 
better  agreement  with  the  data.  A  review  of  reference  (7)  also 
reveals  that  some  difficulty  was  encountered  in  reducing  the  measured 
drag  coefficients,  which  were  obtained  at  angles  of  attack,  to  equiva¬ 
lent  zero  angle-of -attack  coefficients.  It  should  be  noticed,  how¬ 
ever,  that  good  agreement  was  obtained  over  nearly  the  same  Reynolds 
number  range  between  calculated  and  measured  drag  coefficients  pre¬ 
sented  in  figure  2.  The  data  in  figure  2  will  be  discussed  in  more 
detail  later. 

The  two  groups  of  data  in  figure  1,  which  were  reported  in 
reference  (8),  are  at  Reynolds  numbers  near  the  free-stream  transi¬ 
tion  Reynolds  number.  The  value  of  the  free-stream  transition 
Reynolds  number,  as  obtained  from  reference  (8),  is  indicated  by 
the  side  of  each  curve.  These  transition  Reynolds  numbers  were  also 
used  in  the  calculations.  All  data  and  calculations  at  Reynolds 
numbers  below  the  indicated  transition  Reynolds  number  are  for  com¬ 
pletely  laminar  boundary  layer  over  the  body.  Data  shown  at  Reynolds 
numbers  above  transition  Reynolds  number  had  transition  occurring  at 
some  point  on  the  body  and,  hence,  both  laminar  and  turbulent  flow 
were  present.  The  agreement  between  the  calculated  drag  coefficients 
and  the  experimental  data  for  these  higher  Reynolds  numbers  is  good 
for  both  the  Mach  number  9  and  the  !?*■..!;  number  13  cases. 

The  comparison  between  calcv"  ~  J  drag  coefficients  and  data 
presented  in  reference  (9)  is  shov,..  in  figure  2.  These  data  cover 
a  Reynolds  number  range  of  approximately  7  x  10^  to  5  x  106  and  were 
obtained  at  a  free-stream  Mach  number  of  approximately  12.8.  The 
model  used  for  this  test  was  a  slightly  blunted,  8-degree  half -angle 
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cone  with  a  nose-to-base  radius  ratio  of  0.035.  Calculations  at  a 
Reynolds  number  of  approximately  7  x  10^  indicate  that  approximately 
11  percent  of  the  total  drag  is  a  result  of  viscous  interaction 
effects.  Viscous  interaction  effects  which  are  included  in  the 
calculations  represent  increases  in  the  pressure  drag  component  and 
the  friction  drag  component  due  to  induced  pressures  and  increases 
in  the  friction  drag  component  due  to  transverse  curvature.  The 
local  calculated  viscous  interaction  parameter  Ye  at  the  back  end 
of  the  cone  has  a  value  of  2.03  with  an  average  value  of  the  viscous 
interaction  parameter  over  the  frustum  portion  of  the  cone  of  2.5. 


Hie  highest  Mach  number  for  which  the  calculations  were  compared 
with  data  is  approximately  15.  This  comparison  is  shown  in  figure  3 
Hie  data  shown  in  this  figure  were  taken  from  reference  (8).  The 
Reynolds  number  range  for  these  data  extends  from  0.7  to  13  x  10°, 


At  free-stream  Reynolds  numbers  based  on  model  length  greater  than 
5  x  10^,  it  was  assumed  that  transition  occurred  on  the  model.  For 


this  configuration  and  Mach  number,  a  transition  Reynolds  number  of 
5  x  106  was  reported  in  reference  (8).  Again,  over  the  entjre 


Reynolds  number  variation  which  considered  for  some  cases  all  laminar 


flow  on  the  body,  while  for  other  conditions  it  considered  both 


laminar  and  turbulent  flow  existing  on  the  body,  it  is  seen  that  the 
agreement  between  the  data  and  calculations  is  good.  Figure  4  is 
presented  to  show  the  applicability  of  the  program  for  calculating 
drag  coefficients  for  cones  over  a  wide  range  of  nose-to-base  radius 
ratios.  The  data  shown  in  figure  3  were  obtained  from  reference  (10). 
Hie  configuration  used  to  obtain  these  data  was  a  10-degree  half¬ 
angle  cone,  with  nose-to-baso  radius  ratio  varying  from  zero  to 
approximately  0.76.  Again,  the  agreement  between  the  calculated  drag 
coefficients  and  the  measured  values  is  seen  to  be  good. 


LOCAL  FLOW  PROPERTIES 

It  is  often  of  interest  to  have  the  longitudinal  distribution 
of  local  flow  properties  which  exist  at  the  outer  edge  of  the  boundary 
layer.  Also  of  interest  is  the  distribution  of  the  integrated 
boundary-layer  characteristics,  such  as  the  momentum  and  displacement 
thickness.  It  was  mentioned  that  values  for  these  quantities  are 
calculated  along  the  length  of  the  cone.  As  discussed  by  Wilson 
(refs.  (5)  and  (6)),  slightly  blunting  a  slender  cone  flying  at 
hypersonic  speeds  can  significantly  affect  the  values  of  the  local 
flow  properties  over  that  for  a  sharp  cone.  Over  the  portion  of 
the  cone  experiencing  a  strong  entropy  gradient,  the  local  Mach 
number  will  be  much  below  the  value  for  an  invlscid  sharp  cone.  The 
local  Mach  number  can  be  used  to  indicate  the  "swallowing  length," 
or  the  distance  required  for  the  boundary  layer  to  "swallow"  the 
variable  entropy  layer.  The  swallowing  length  can  be  defined  as 
that  distance  required  for  the  local  Mach  number  to  reach  some  arbi¬ 
trarily  chosen  fraction  of  the  local  inviscid  cone  Mach  number,  say 
99  percent.  To  illustrate  the  effect  nose  blunting  has  on  the  local 
Mach  number,  the  distribution  of  the  local  Mach  number  at  the  outer 
edge  of  the  boundary  layer  on  a  5-degree  half-angle  cone  has  been 
calculated  and  presented  in  figure  5.  These  calculations  were 
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performed  for  a  sharp  cone  and  for  two  slightly  blunted  cones.  Hie 
free-stream  Mach  number  and  slant-length  Reynolds  number  used  in 
these  calculations  were  12  and  1  million,  respectively. 

To  illustrate  the  effect  of  blunting  on  integrated  boundary- 
layer  characteristics,  the  results  of  calculations  of  the  local 
momentum-thickness  Reynolds  number  distribution  is  shown  in  figure  6. 

Prom  these  two  figures  it  is  seen  that  slight  blunting  reduced 
the  value  of  both  the  local  Mach  number  and  momentum-thickness 
Reynolds  number  at  any  given  point  along  the  surface  of  the  cone. 
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FIG  .3  VARIATION  OF  DRAG  COEFFICIENT  WITH  REYNOLDS  NUMBER  M  =15,  0  =  9  DEGREES 
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APPENDIX  A 

This  appendix  includes  a  description  of  the  details  of  the 
program  as  it  is  presently  coded  for  use  on  an  I3M  7090*  Included 
in  this  description  is  the  required  input  and  the  input  format  for 
the  program.  Also  included  is  a  cross  listing  of  the  symbols  used 
in  the  text,  with  those  appearing  in  the  Fortran  listing  and  those 
appearing  in  the  program  output.  Further,  there  is  a  listing  of 
the  definitions  of  the  symbols  used  in  the  program  output. 

This  program  is  coded  in  Fortran  IV  and  nans  under  the  IBSYS  - 
IBJOB  monitor  version  13.  Pertinent  comment  cards  have  been  inserted 
into  the  program  listing  to  aid  in  relating  program  statements  to 
equation  in  the  text  or  given  references.  Tnis  was  done  to  make  any 
desired  modification  of  the  program  as  easy  as  possible. 

This  program  will  compile  and  run  the  test  case  (output  in 
Appendix  B)  on  an  IBM  7090  in  1.53  minutes. 
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Card 

IDEN 

2A6 

Card 

RNBB 

Sia.6 


SYMBOLS 

IDEN 

SIQMA 

EMIN 

TIN 

TW 

RORPIN 


RNRB 

RB 

RLTR 


Runs  may 


PROGRAM  INPUT 


1. 


SIGMA 

EMIN 

TIN 

TW 

RORPIN 

E12.6 

E12.6 

El  2. 6 

E12.6 

E12.6 

2. 

RB 

RLTR 

E12.6 

E12.0 

Definition  of  Symbols 

For  Identification  purposes  only 
Cone  half -angle  (degrees) 

Free -stream  Mach  number 

Free-stream  temperature  (degrees  Rankine) 

Wall  temperature  (degrees  Rankine) 

Free-stream  Reynolds  number  based  on  axial  length  (inches) 
or  free-stream  pressure  (atmospheres) 

Ratio  of  nose-to-base  radii 
Base  radius  (inch) 

Free-stream  transition  Reynolds  number  based  on  axial 
length  (inch) 

be  stacked  with  a  blank  card  after  last  run. 
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INTERNAL  AND  EXTERNAL  COMPUTER  SYMBOLS 


Computer 

Symbol 

Symbol 

Me 

EMC 

Meo 

EMIN 

Y 

GAM 

Pc 

PC 

PIN 

rB 

RB 

RLro(tr) 

RLTR 

*nAb 

RNRB 

a 

SIGMA 

T 

■‘■CO 

TIN 

Tw 

TW 

Tc 

TC 

°f 

CF 

% 

CIBAR1 

6 

DELTA 

Mi 

Ml 

R*c 

RC 

Rsl 

RES 

w 

CD 

o 

RTHC 

R@1 

RTH1 

S 

S 

*1 

T1 

X 

X 

% 

CDB 

cx 

CDFL 

cn 

CDFT 

Df  T 

^Df( total ) 

CDFTTL 

CDN 

S, 

CDPL 

CDPT 

cDp( total ) 

CDPTTL 

Cj)(  total ) 

CDTTL 

CD( total) 

CDTTLB 

Output 

Symbol 

MC 

M(  INF) 

PC 

P(  INF) 

RB 

RL(TR) 

RN/RB 

SIGMA 

T(INF) 

T(WALL) 

TC 

CF 

CIBAR1 

DELTA 

Ml 

RXC 

RSI 

RTHC 

RTH1 

S 

T1 

X 

CD( BASE) (VACUUM) 

CDF(L) 

CDF(T) 

CDF/jOTAL) 

CD(NOSE) 

CDP(L) 

CDP(T) 

CDP(  TOTAL) 

CD( TOTAL)  excluding  base  drag 
CD( TOTAL)  including  base  drag 
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Internal  and  External  Computer  Symbols  -  continued 


Computer 

Output 

Symbol 

Symbol 

Symbol 

ACDf(IP) 

DCDPIP 

D  CDP(IP)(L) 

ACoplTC) 

DCDPTC 

D  CDF(TC )(L) 

ACDp(IP) 

DCDPIP 

D  CDP(IP)(L) 

*A 

CIBARA 

CHIBAR  AVG(L) 

INTERNAL 

Computer 

COMPUTER  SYMBOLS 

Symbol 

Symbol 

6/e 

DOVTH 

Defined  EQ.  (12)(  Ref .( 6. ) ) 

d(6/e-H)/dM-, 

DTHDM1 

Defined  EQ.  (13)(  Ref .( 6. ) ) 

Y-l 

GM1 

Y+l 

GP1 

IU 

OMEGA 

Shock  wave  angle  for  sharp  cone 

wa 

OMEGAS 

Local  shock  wave  angle  for  blunt  cone 

Pc/P„ 

PC  PIN 

Ratio  cone  pressure  to  free-stream 
pressure 

rr 

PI 

rBC 

RBC 

Reynolds  number  based  on  cone  properties 
and  tip  radius  of  curvature 

rn 

RN 

Nose  radius 

Sf£ 

ST 

Distance  along  the  surface  of  nose 
from  stagnation  to  tangency  point 

T'Ai 

TPT1 

Defined  EQ.  (4X  Ref .(  5- ) ) 

TT1 

Stagnation  temperature  at  the  outer 
edge  of  boundary  layer 

U'/htu 

XMU1MW 

u'/ui 

XMPM1 

Ul/'Mc 

YMT  MQ 

*«*•*  »v 

F1 

FI 

Function  defined  by  EQ,  ( 1 3 ) ( Ref . ( 5 . ) ) 

§ 

XRRR 

RcRec/RB^  “  X0/*§ 

2f/F2 

YRRR 

Integrand  in  EQ.  (20) (Ref .( 5. ) ) 

f2 

F2 

Function  defined  by  EQ.  (l6)(Ref .( 5* ) ) 

SIN^ujs 

SINWS2 

C0S2iua 

C0SWS2 

TAN2uja 

TANWS2 

TAN2uj 

TANW2 
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SYMBOLS 

P(INP) 
T(  INF) 
M(INF) 

PC 

TC 

MC 

SIGMA 
T(WALL  ) 
RN/RB 

RB 

RL(INF) 

RL(TR) 


D  CDP(  IP )  ( L ) 


D  CDF(  TC  )  ( L ) 


D  CDF(  IP )  ( L ) 


CHIBAR  AVG(L) 


CDF(L) 


CDP(L) 


CD(NOSE) 
CDF (TOTAL) 


CDP( TOTAL) 


PROGRAM  OUTPUT 

Definition  of  Symbols 

Free-stream  pressure  (atmospheres) 

Free-stream  temperature  (degrees  Rankine) 

Free-stream  Mach  number 

Inviscid  sharp  cone  pressure  (atmospheres) 

Inviscid  sharp  cone  temperature  (degrees  Rankine) 
Inviscid  sharp  cone  Mach  number 

Cone  half -angle  (degrees) 

Wall  temperature  (degrees  Rankine) 

Ratio  of  nose-to-base  radii 

Base  radius  (inch) 

Free-stream  Reynolds  number  based  on  axial  length 
Free-stream  transition  Reynolds  number  based  on 
axial  distance  from  stagnation  point 

Incremental  increase  in  the  pressure  drag  coefficient 
due  to  induced  pressure  for  the  laminar  portion  of 
the  boundary  layer 

Incremental  Increase  in  the  skin-friction  drag 
coefficient  due  to  transverse  curvature  for  laminar 
portion  of  the  boundary  layer 

Incremental  increase  in  the  skin-friction  drag 
coefficient  due  to  Induced  pressure  for  laminar 
portion  of  the  boundary  layer 

Average  value  of  viscous  Interaction  parameter, 

(X)  for  laminar  portion 

Skin-friction  drag  coefficient  for  laminar  portion 
of  boundary  layer 

Pressure  drag  coefficient  for  area  wetted  by 
laminar  boundary  layer 

Pressure  drag  coefficient  for  spherical  nose  cap 

Total  skin-friction  drag  coefficient  (sum  of  all 
skin-friction  drag  coefficient  components) 

Total  pressure  drag  coefficient  (sum  of  all 
pressure  drag  coefficient  components) 


CD( BASE )( VACUUM )Base  drag  coefficient  assuming  vacuum  at  the  base 
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CD( TOTAL ) 
excluding 
base  drag 

CD( TOTAL) 
including 
base  drag 


X 

S 

Ml 

CP 


Sum  of  all  drag  coefficient  components  except 
base  drag  coefficient 


Sum  of  all  drag  coefficient  components 


Distance  along  cone  surface  measured  from 
hypothetical  apex  (inch) 

Distance  along  cone  surface  measured  from 
stagnation  point  (inch) 

Local  Mach  number  at  outer  edge  of  boundary  layer 
Local  skin -friction  coefficient 


T1 

DELTA( 6 ) 
RXC 

RTHC 

RTH1 

RSI 

CHIBARI 


Local  temperature  at  outer  edge  of  boundary  layer 
(degrees  Rankine) 

Boundary -layer  thickness  (inch) 

Reynolds  number  based  on  inviscid  sharp  cone 
properties  and  X 

Reynolds  number  based  on  inviscid  sharp  cone 
properties  and  boundary-layer  momentum  thickness 

Reynolds  number  based  on  local  properties  and 
boundary-layer  momentum  thickness 

Reynolds  number  based  on  local  properties  at 
outer  edge  of  boundary  layer  and  S 

Local  value  of  the  viscous  interaction  parameter  x 
defined  by  _  m5  r— 


VrsT 


where  C  «  pw 

HI  / 


Tw 

*1 


<3 
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PROGRAM  LISTING 

SIBFTC  DRAG  DD.DECK 

THIS  PROGRAM  IS  CODED  IN  FORTRAN  IV  AND  RUNS  UNDER  THE  IBSYS- 
IBJOB  MONITOR  VERSION  13. 

MAIN  PROGRAM 

INPUT  CONSIST  OF  2  CARDS  PER  ROUND. 

FIRST  CARD  -  FORMAT ( 2 A6 , 5E1 2 .6 ) 

IDEN  -  IDENTIFICATION  (MAX.  12  CHARACTERS) 

SIGMA  -  CONE  HALF-ANGLE  (DEGREES) 

EMIN  -  FREE  STREAM  MACH  NO. 

TIN  -  FREE  STREAM  TEMPERATURE  (RANKINE) 

TW  -  WALL  TEMPERATURE  (RANKINE) 

RORPIN-  FRFE  STREAM  REYNOLDS  NUMbER  BAsED  ON  AXIAL  LENGIH  IN 
INCHES  (OR)  FREE  STREAM  PRESSURE  (ATMOSPHERES) 

NOTE  -  IF  RORPIN  LESS  THAN  100.  PROGRAMS  ASSUMES  PRESSURE. 

IF  RORPIN  GREATER  THAN  100.  PROGRAM  ASSUMES  REYNOLDS  NUM. 

SECOND  CARD  -  FORMAT ( 6E12. 6 ) 

RNRB  -  RATIO  NOSE  TO  BASE  RADIUS. 

RB  -  BASE  RADIUS  (INCHES) 

RLTR  -  FREE  STREAM  TRANSITION  REYNOLDS  NO.  bASED  ON  AXIAL  LENGTH 

IF  AN  ALL  TURBULENT  RUN  IS  DESIRED  SET  RLTR  EQUAL  TO  A  SMALL 
VALUE.  SAY  100. 

ROUNDS  MAY  BE  STACKED  -  LAST  ROUND  SHOULD  BE  FOLLOWED  WITH  A  BLANK  CARD. 

PROGRAM  OPERATING  MODES 
M0DE=1  -  LAMINAR 
M0DE=2  -  LAMINAR-TURBULENT 
M0DE=3  -  TURBULENT 


PROGRAM  CONTROL  VARIABLES 

COMMON  DELM.DELRC.DELW, IQUI T , I SK I P , I TRbLE , ITRY , K . KODE .KEY » 

*  LINES, MODE, N.Nl.XLIMIT 
PROGRAM  INPUT  VALUES 

COMMON  EMC.EMIN.GaM, IDENt  2) .PC.PIN.Rb.RLIR ,RNRB , S IGMA, T IN »TW 
PROGRAM  NO  ND I  ME  N  S I  ON  E  D  OUTPUT  VARIAbLES 
COMMON  CC3,CDF,CDFL,CDFT,Ci/Fl  I  L ,  COw ,CDP ,  CDPL  »CDP  I  ,CDPI  IL.CDT  It, 
*CDTTLB»DCDFIP»DCDFTC>DCDPIP»CIBARA,TC 
PROGRAM  DIMENSIONED  0U1PUI  VARIABLES 
COMMON  CF ( 500 ) »CI BARI ( 500 ) ,DELTA(500) ,EM1(500) »RC(500) » 

*RES ( 500 ) ,RTHC(500) ,RTH1 (500) ,S(50C)  , T 1 ( 500 ) »X(500) 

PROGRAM  NONDIMENSIONED  INTERNAL  VARIABLES 
COMMON  CONST, CPC, C0SWS2 ,DELS,DIN ,D0VTH,DTHDM1 ,GM1 ,GP1, 

*  OMEGA, PCPIN, PI ,  RBC,RCTR,RN,SI NWS2 ,SMAX  »SMAXT  » 


" ST  ,  T  ANW2 , T  AN  W32,TPTl,T71»xri,XF2»  XMAX , XMPM1 , XMU1MW , XM1MC 


PROGRAM  DIMENSIONED  INTERNAL  VARIABLES 


COMMON  FI ( 500 ) .OMEGAS ( 500 ), XRRR (  750)  ,YRRR(75C) 


SINCE  A  CONSTANT  PRESSURE  ON  I  HE  CONICAL  SURFACE  EQUAL  TO  THE 
INVISCID  SHARP  CONE  PRESSURE  IS  ASSUMED,  WE  SET  P1=PC  AND  Pl/PIN 
=PC/P I N  USING  AN  EQUIVALENCE  STATEMENT. 

EQUIVALENCE  (PC, PI ) , ( PCPI N ,P1 PIN ) 
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READ  INPUT  FIRST  CARD 
5  READ ( 5  » 1 0 )  IDEN, SIGMA, EMIN'  iIN.TW.RORPIN 
10  FORMAT (2A6.5E12.6) 

IFISIGMA.EQ.O. )  STOP 

READ  INPUT  SECOND  CARD 
READ! 5.15)  RNRB  »RB » RLTR 
15  FORMAT ( 6E12  «  6 ) 

SETUP  PROGRAM  CONTROLS.CONSTANTS  AND  UNIT  CONVERSIONS 
I TRBLE=0 
KEY  =  0 

RN=RMR3*R6 

SET  GAMMA  =  1.4 
GAM= 1.4 
GM1=GAM-1 . 

GP1=GAM+1 • 

PI=3. 14159 

SIGMA=SIGMA/57. 29578 
DELW=4./57. 29578 
IQU I  T  =  0 

CALCULATE  GEOMETRIC  CHARACTERISTICS  OF  SPHERICALLY  BLUNTED  CONE 
SMAX=(RN*(PI/2.-SIGMA)+(RB-RN*SIN(PI/2.-SIGMA) ) /SIN (SIGMA) ) 
ST=RN*(PI/2. -SIGMA) 

AXIAL-RN*( 1. -SIN (SIGMA) ) + ( Rb-RN*SIN ( P I /2.-S IGMA) ) *COS( SIGMA ) 

1 /S  IN ( S I GMA ) 

P I N=RORP I N 

IFIRORPIN.LT. 100. )  GO  TO  25 

PIN=RORPIN/(222.*EMIN*AXIAL*l.E6*(TIN+198.6)/T IN**2) 

CALCULATE  INVISCID  SHARP  CONE  PROPERTIES 
25  CALL  PCMCTC 

TANW2=TAN(OMEGA)**2 
MODE= 1 

IF (RLTR.GT.O. )  MODE=2 
GO  TO  (40.30) .MODE 

CALCULATE  S  VALUE  WHERE  TRANSITION  WILL  OCCUR  IF  RLTR  NOT  ZERO 
30  SMAXT=RN*(PI/2.-SIGMA)+(RLTR*TIN**2*l.E-6/(222.*EMIN*PIN*( I IN+ 
1198.6) )-RN*( 1. -SIN (SI GMA) ) ) /COS (SIGMA) 

33  IF ( SMAXT.LT .ST )  SMAXT=ST 
IFISMAXT.LT. SMAX)  GO  TO  35 
MODE= 1 
GO  TO  40 
35  SMAX=SMAXT 

LAMINAR  BOUNDARY  LAYER  CALCULA1IONS 
40  CALL  LAM  I  NR 

IF( ITRBLE.GT.O)  GO  TO  5 

SETUP  PROGRAM  CONTROLS  AND  CONSTANTS 
DELS=S ( 1 ) -ST 
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Nl  =  l 

C0FTTL=0. 

CDPTTL=0. 

IF (M0DE.EQ.2.AND.N.EQ.1 )  GO  TO  60 

CALCULATE  VISCOUS  INTERACTION  tFFECTS  ANU  DRAG  COEFFICIENT 
COMPONENTS  FOR  LAMINAR  PORTION 

CONST=RB**2*(RNRB**2*(2.*( 1 .-S IN ( SIGMA) ) -COS ( S I GMA ) **2 /S I N ( SIGMA) ) 
1+1 » /SI N ( S 1 GMA ) ) 

XF1*. 8992-. 1049*EMIN#SI GMA 

XF2=.6192-,3822»EMIN*SIGMA+.0798#EMIN*#2*SIGMA*#2 
DIN=.968#TW/(EMIN**2*TIN)+.145*GM1 
CALL  SUMUP 
C 

C  SUM  LAMINAR  DRAG  COMPONENTS.  SETUP  PROGRAM  CONTROLS  AND  ETC. 

Nl-N 

CDFL=CDF 

CDPL«=CDP 

CDFTTL=CDFL+DCDFIP+UCOFTC 
CDPTTL=CDPL+DCDPIP 
IF (MODE.EQ. 1 )  GO  TO  70 
N1«N 

60  SMAX=RN*(PI/2.-SIGMA)+(RB-RN#SIN<PI/2.-SIGMA))/SIN(SIGMA) 
IF(N.EQ.l)  MODE=3 

XMAX=SMAX-RN# <  P I / 2 . -S I GMA } +RN / TAN ( S I GMA ) 
RFINAL=222.*EMC#PC#(TC+198.6)#XMAX*1.E6/TC**2 
DELRC* (KFINAL-RC ( N ) )*XMAX/< (XMAX-XtN) )*100.) 

C 

C  TURBULENT  BOUNDARY  LAYER  CALCULATIONS 
CALL  TURB 

IF (  ITRBLE.GT.O)  GO  TO  5 
C 

C  CALCULATE  DRAG  COEFFICIENT  COMPONENTS  FOR  TURBULENT  PORTION 
DELS=S(N1+1)-S(N1 ) 

CALL  SUMUP 

CDFT=CDF 

CDPT=CDP 

CDFTTL=CDFTTL+CDFT 

CDPTTL=CDPTTL+CDPT 

C 

C  CALCULATE  DRAG  COEFFICIENT  FOR  SPHERICAL  NOSE  CAP 

70  CDN=  (  (6.*EMIN**2/5.  ) **3 •  5* (  6*  /  ( 7.*£MIN*"*2-1 . )  ) **2 •  5—1.  )*RNRB#*2* 

1 ( 1. -SI N ( SIGMA) **4)/(1.4*EMIN**2) 

C 

C  CALCULATE  BASE  DRAG  COEFFICIENT  (ASSUME  VACCUUM) 

CDB=1./(GAM/2.*EMIN**2) 

C 

C  SUM  UP  LAMINAR  AND  TURBULENT  COMPONENTS  TO  FORM  TOTALS 
CDTTL=CDFTTL+CDPTTL+CDN 
CDT  TLB=CDTTL+C0t3 
C 

C  PRINT  RESULTS 

CALL  OUTPUT 
GO  TO  5 
END 
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IBFTC  CONPRO  00 

THIS  SUBROUTINE  CALCULATES  INVISCID  SHARP  CONE  PROPERTIES 
SUBROUTINE  PCMCTC 

COMMON  DELM.DELRC.DELW,  IQUIT,  I  SKIP,  ITRbLt,  I  TRY »K , KODE .KEY » 

*  L I NES  »M0DE  .N.Nl.XLIMIT 

COMMON  EMC  >  EMIN  .GAM  ,  I OEN  (  2  )  .PC  ,P  I N ,  Ru  ,RL  I  R  ,RNRB  ,b  IGMA,  I  IN»IW 
COMMON  CDb»COF  »CDf-  L  ,CDFT  »CUFT  TL  »CUN  >C0r  ,CDPL ,CDPT .CDPTTL .CDTTL , 
*CDTTLB»DCDFIP»DCL'FTC»DCDPIP»CIBARA»TC 
COMMON  CF( 500) »CI BARI (500) .DELTA (500) .EMI ( 500) »RC(500) . 

*RES<500) ,RTHC(500) .RTH1 !5005 .S(500) .Tl (500) »X(500) 

COMMON  CONST. CPC >COSWb2, DELS, 0  I N.DOVTH.DTHDMl.GMl.GPl, 

*  OMEGA, PCPIN, PI  ,  RBC.RCTR.RN.SI NWS 2 .SMAX.SMAXT , 

*ST  »TANW2 .TANWS2 , i PT 1 , TT 1 , XF1 , XF2 ,XMAX , XMPM1 .XMU1MW , XM1MC 
COMMON  FKSOOi  ,OMEGAS(  500)  .XRRR(75o)  »YRRR(75o) 

EQUIVALENCE  (PC,P1),(PCPIN,P1PIN) 

XKC-EM I N*S IGMA 
GP3=GAM+3. 

XKS=GP1*XKC/GP3+SQRT( ( GP 1/GP3 ) **2*XKC**2+2 ./GP3 ) 

CALCULATE  CONE  TEMPERATURE  FROM  REF.  3 
I F ( EMI N*SIN ( SIGMA ) «GT ♦ 1 . )  GO  TO  10 
TCTIN=l.+.35*( EMI N*SI N( SIGMA ) )**1.5 
GO  TO  25 

10  TCTIN=( l.+EXP(-l.-1.52*EMI N#S I N ( S I GMA ) ) )*( 1 .  +  ( EMI N*S1 N ( S IGMA )  1**2/ 
14.  ) 

CALCULATE  CONE  PRESSURE  FROM  REF.  4 
25  PCPIN=l.  +  .7*EMIN**2*4.*SlN(S>IGMA)**2*(2.5+8.*bQRl  (  EMIN**2-1 »  )  *bIN 
1  (SIGMA)  ) /(1.+ 16. *SQRT < EMI N**2-l.)*SlN( SIGMA)  ) 

PC=PCPIN*PIN 

CALCULATE  CONE  MACH  NUMBER  FROM  REF.  3 
UCUIN=COS ( SI GMA )*SQRT(1. -SIN (SIGMA) /EMIN) 

EMC=UCUI N*EM IN/SQRT ( TCT IN ) 

40  TC  =  TCT I N*TI N 

CPC= (PCP IN-1. )/ (GAM/2. *EMIN**2> 

CALCULATE  SHOCK  WAVE  ANGLE  FOR  SHARP  CONE 
OMEGA=XKS/EMI N 
RETURN 
END 
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SIBFTC  LAMNR.  DD.DECK 
SUBROUTINE  LAMINR 

C  THIS  ROUTINE  CALCULATES  LOCAL  PROPERTIES  FOR  THE  LAMINAR  PORTION 
COMMON  DELM.DELRC.DELW, IQU I T , I SK I P , I TRbLE , I  TRY ,K , KODE ,KEY * 

*  LINES, MODE, N.Nl.XLIMIT 

COMMON  EMC,  EMIN, GAM,  I  DENI  2)  ,  PC  ,P  I N ,  Rb  ,RLTR  ,RNRB  ,S  IGMA,  T I N ,  TW 
COMMON  CDb,CDF,CUFL,CUFT,CuFTTL,CUN,CDP,CDPL,CDPT ,CDPTTL,CDTTL, 
#CDTTLB,DCDFIP,DCDFTC,OCDPIP,CIbARA,TC 
COMMON  CF( 500) ,CIBAR1 ( 500) .DELTA (500) , EMI ( 500 ) , RC ( 500 > » 

*RES(500)  »7THC(500) »RTH1(500) , S ( 500 )  , T 1 ( 500 ) , X ( 500 ) 

COMMON  CONST, CPC, COSWS2, DELS, D I N,D0VTH,DTHDM1,GM1,GP1, 

*  OMEGA, PCPIN, PI  ,  RBC,RCTR,RN,SIiV;.'S2»SMAX,SMAXT, 

*ST ,TANW2 ,TANWS2 ,TP  T1 »  TT1 ,XF1 , XF2 ,XMAX,XMPM1 ,  XMU1MW , XM1MC 

COMMON  FI ( 500  )  .OMEGAS (500),XRRR(750)  ,YRRR(750) 

EQUIVALENCE  ( PC  ,P 1 ) , ( PCPI N  ,  PI P  IN ) 

FUNCTION  STATEMENT  FOR  LINEAR  INTtRPOLAT I  ON 
XINT(XPT,X1»X2»Y1 , Y2 ) =Y1+ ( Y2-Y1 ) * ( XPT-X1 ) / ( X2-X1 > 

I FEW=0 
I MANY=0 

IFIRN.GT.O. )  GO  TO  1 

LAMINAR  BOUNDARY  LAYER  CALCULATIONS  FOR  SHARP  NOSED  CONE 
CALCULATE  100  POINTS  OR  EVERY  .01  INCHES  ALONG  BODY 
195  DELX=SMAX/100. 

I F ( DELX.LT . .01 )  DELX=.01 

X ( 1 ) =0 . 

N=  10 1 

CALCULATE  T ( PR IME ) / T ( C ) 

TPTC=1.+.076*GM1*EMC#*2+.481*(TW/TC-1. ) 

CALCULATE  MU ( PRI ME ) /MU ( C ) 

XMPMC=  (TC+198.6)  *TPTC*SQRT  (TPTC)/(TPTC*TC+198.6) 

CALCULATE  F(2)  EQ.  (23),  REF.  5 
F2=.441*XMPMC/TPTC 
DO  200  J=l,101 

CALCULATE  REYNOLDS  NUMBER  bASED  ON  CONE  PROPERTIES  AND  X 
RC(J)=222.#EMC*PC#X(J)*(TC+198.6)*1.E6/TC*#2 

CALCULATE  R  THETA ( C )  EQ.  (24),  REF  5 
RTHC ( J ) =SQRT (F2*RC(J)/3.) 

CALCULATE  LOCAL  SKIN  FRICTION  COEFFICIENT  EQ.  (26),  REF  5 
CF ( J  )  =F2/RTHC ( J ) 

CALCULATF  BOUNDARY  LAYER  THICKNESS  EQ.  (5),  REF  5 
DELTA(J)=(6.1*TPTC+2.9)*RTHC(J)#X(J)/RC(J) 

EMI ( J) =EMC 
T1 ( J ) =TC 
RTH1 ( J ) =RTHC ( J ) 

RES ( J ) =RC ( J ) 

CE=(T1(  J)  +  198.6)*SQRT(Tw/Tl(J)  )  /  (  Tv<+198.6  ) 
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C  CALCULATE  LOCAL  CHI  BAR  PARAMETER 

C I  BARI ( J ) = EMC* *3* SORT ( CE ) /SORT ( RES ( J )  ) 

S< J)=X( J) 

I F  t  X ( J) .GE.SMAX)  GO  TO  150 
X< J+1)=X( Jl+DELX 
200  CONTINUE 
J=10l 
GO  TO  150 
C 

C  LAMINAR  BOUNDARY  LAYER  CALCULATIONS  FOR  BLUNTED  CONE. 

1  I TRY=0 
XL  I  MI T  =  1 . 

OMEGAS! 1 )=PI /2» 

CALCULATE  REYNOLDS  NUMBER  BASED  ON  CUNE  PROPERTIES  AND  NOSE  RADIUS 
RBC=222.*PC*EMC*RN*( ( TC+198 . 6 ) /TC»*2 » * ' . t6 

CALCULATE  REYNOLDS  NUMBER  BASED  ON  CONE  PROPERTIES  AND  BOUNDARY- 
LAYER  MOMENTUM  THICKNESS  AT  POINT  OF  TAN6ENCY . 
RTHCN=6.48*SQRT(1.5*RBC) 

CALCULATE  REYNOLDS  NUMBER  oASED  ON  CONE  PROPERTIES  AND  DISTANCE 
ALONG  CONE  SURFACE  MtASUREU  FRUi>.  mPEX  10  I  Ai^GENCY  POINT. 
RCN=222.*PC*EMC*RN*1.E6*( TC+198.6 )/( TC**2*TAN ( SIGMA ) ) 

CALCULATE  BOUNDARY-LAYER  GROWTH  XI (NOSE)  OVER  THE  NOSE  BY  METHOD 
OF  REF.  12  IN  REF.  5. 

XRRRN=43 . /SORT ( RBC ) 

5  I QU I T=0 
K0DE=0 

DO  100  1=1,440 

START  1 1  G  WITH  PI/2  VARY  OMEGA ( S )  IN  INCREMENTS  OF  DELTA  OMEGA(S) 

10  OMEGAS!  I  i-'  =OMEGAS (  I  )  -DELW 
20  SINWS2=S1  , OMEGAS ( I +1 ) ) **2 
TANWS2=TAN ( OMEGAS ( I +1 ) ) **2 
COSWS2=COS ( OMEGAS ( I +1 ) )**2 

CALCULATE  LOCAL  MACH  NO.  FROM  EO.  (11)  REF  5. 

TEMP1=(GP1/ <2.*GAM*EMIN**2*SINwS2-GMl ) )**( l./GAM) 
TEMP2=GP1*EMIN**2*SINWS2*( 2 . +GM1*EMI N**2 > 

TEMP 3=2. *( 2 .+GM1*EMIN**2*SINWS2 ) 

EM1TMP=SQRT( ( ( TEMP 1*T EMP2 /TEMP3 ) /PI P I N** ( GMl /GAM ) -1 . ) / ( GM1/2 . ) ) 

CALCULATE  XI  FROM  EO.  (21),  REF  5. 

TEMP4=SQRT ( ( 2.+GMl*EMIN**2 ) / ( 2 .+GM1*EM1TMP**2 ) ) 

XRRR ( I )  =  .1755*EMIN*TEMP4/ ( PIP  I N*SIN ( S IGMA ) *EM1TMP* ( TANWS2-TANW2 ) ) 
IF ( I .LE.l.OR. (XRRR( I ) -XRRR ( 1—1 ) ) .Lt.XLlMIT)  GO  TO  30 
DELW=PELW/? . 

GO  TO  10 

CALCULATE  F (  1 )  EO.  (13),  REF  5. 

30  TEMP5=2.+GM1*EM1TMP**2-6.41 

TEMP6=EM1TMP*#2*TAhWS2/ ( 2.*COSwS2* ( TAN*S2-TANw2 ) ) 
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TEMP7=(2.*GAM*EMIN**2*SINWS2-GM1 )*( 2.+GMl*EMIN**2*SINWS2 ) / 

1( EMIN**2*SINWS2-1. )**2 
TEMP8=1.+GM1*EM1TMP**2 
F1TMP=TEMP5/(TEMP6*TEMP7-TEMP8  ) 

C 

C  CALCULATE  T ( 1 )  EQ.  (18).  REF  5. 

T1 TMP=TC* ( 2 .+GM1*EMC**2 ) / < 2 .+GM1*EM1TMP**2 ) 

CALCULATE  T ( PR  I  ME ) /T ( 1 )  EQ.  (A),  REF  5. 
TPT1=1.+.076*GM1*EM1TMP**2+.481*(TW/T1TMP-1. ) 

CALCULATE  MU( PRIME) /MU( 1 )  FROM  SUTHERLANDS  VISCOSITY  LAW. 
XMPM1=(T1TMP+198.6) *TPT1** 1 . 5/ < TPT1*T 1 TMP+198 .6 ) 

CALCULATE  MU(1)/MU(C)  FROM  SUTHERLANDS  VISCOSITY  LAW. 
XMlMC=(TC+198.6)*(TlTMP/TC)**i»5/(TlTMP+198«6  ) 

CALCULATE  F ( 2 )  USING  EOS.  (16)  AND  (17).  REF  5. 

F2=.441*EMC*SQRT(  T 1  TMP/ TC  )  *XM1MC*XMP,-H/ (  TPT 1  *£M1TMP*  ( 1.  +  F1TMP)  ) 

CALCULATE  INTEGRAND  IN  Eg.  (20).  REF  5. 

YRRR ( I )=2.*XRRR( I )/F2 
IF(KODE.GT.O)  GO  TO  50 

IF  XI  IS  BELOW  THE  LOWER  LIMIT  (XI (NOSE)),  KEEP  STEPPING  OMEGA(S) 
UNTIL  XI  EQUALS  OR  EXCEEDS  LOWER  LIMIT,  INTERPOLATE  IF  NECESSARY. 
IF(XRRR( I ) .LE.XRRRN)  GO  TO  100 
I  F ( I.GT.l)  GO  TO  40 
DELW=DELW/2. 

GO  TO  10 

40  YRRR  (  I  -1 )  =  X  I  NT  ( XRRRN » XRRR  (  1-1 )  ,XRRR( I  )  ,  YRRR  (  1-1 )  ,YRRR(  I  )  ) 

XRRR ( 1-1 ) =XRRRN 
J  =  0 
KODE= 1 
AREA=0 
50  J=J+1 

USING  TRAPEZOIDAL  RULE,  INTEGRATE  THE  INTEGRAL  IN  EQ.  (20).  REF  5. 
AREA  =  AREA+ ( YRRR ( I ) +YRRR ( I - 1 ) )/2.*(XRRR( I )-XRRR( 1-1 > ) 

CALCULATE  R(C)  EO.  (20),  REF  5 

RC( J)=R6C*(3.*RBC*AREA+l./TAN(SI&MA)#*3)**(l./3. ) 

CALCULATE  X  FROM  R ( C  >  AND  CONE  PROPERTIES. 

X( J)=RC{ J)/(222.*PC*EMC#1.E6*(TC+198.6)/TC**2) 

CALCULATE  S  USING  X  AND  bODY  GEOMETRY. 

S ( J ) =RN* ( P I /2. -S IGMA )  +  ( X ( J )-RN/T  AN ( SIGMA ) ) 

EMI ( J ) =EM1 TMP 

CALCULATE  R  THtTA(C)  USING  R(C),  EQ.  (20)  AND  EQ(21)»  REF  5. 

RTHC ( J ) =XRRR ( I ) *RBC**2/RC ( J ) 

T 1 ( J ) =T1TMP 
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C  CALCULATE  R  THETAtl!  USING  EQ.  (22).  REF.  5. 

RTH1 < J)=RTHC( J)*EM1( J)*(T1 ( J ) +1 98 . 6 ) * ( TC/T 1 ( J ) ) **2/ ( EMC* ( TC+198 . 6 
1  )) 

C 

C  CALCULATE  REYNOLDS  NUMBER  BASED  ON  LOCAL  PROPERTIES  AND  S. 

RES( J)=222.*PI*EM1( J)*S( J)*(T1 ( J)+198.6)*1.E6/T1( J)**2 
C 

C  CALCULATE  LOCAL  SKIN-FRICTION  COEFFICIENT  EQ.  (3).  REF.  5. 
CF(J)=.441*XMPM1/(RTH1(J)*TPT1) 

CE=(T1( J) +198.6 )*SQRT(TW/T1< J) )/ (TW+198.6) 

CALCULATE  LOCAL  CHI  BAR  PARAMETER. 

Cl  BARK J)=EM1( J ) **3*SQRT (CE)/ SORT! RES (J) ) 

Fl( J)=F1TMP 

CALCULATE  BOUNDARY  LAYER  THICKNESS  EQ.  (5),  REF*  5. 

DELTA! J)=(6,1*TPT1+2.9)*RTHC< J)*X< J)/RC( J) 

54  I F ( S ( J ) -SMAX5  100.150. 60 

55  DELW=DELW/5. 

J=J-1 

AREA-: AREA-  (  YRRR  (  I  )  +YRRR  (  1-1 )  )  /2.*(XRRR(  I  )-XRRR(  I  — 1 )  ) 

I QU I T= 1 
GO  TO  10 

CHECKS  TO  CONTROL  NUMBER  OF  LOCAL  VALUES  CALCULATED  OVER  THE 
LAMINAR  PORTION. 

IF  ALL  TUR8ULENT  RUN,  CALCULATE  ONLY  ONE  LAMINAR  POINT. 

IF  LAMINAR-TURBULENT  RUN,  CALCULATE  MINIMUM  15  LAMINAR  POIN'  S 
UNLESS  S ( TR  )  LESS  THAN  SMAX/3. 

IF  ALL  LAMINAR  RUN,  CALCULATE  MINIMUM  25  POINTS  WHEN  FEASIBLE. 

60  IF(MODE.EQ.2.AND. J.EQ.l )  GO  TO  150 
I F ( MODE .EQ.2.AND.J.GT.15)  GO  TO  170 
IFIMODt. EQ.?. AND. S( J) .GT. (SMAX/3. ) }  GO  TO  115 
I F ( MODE .EQ. 1.AND.J.LT.25)  GO  TO  115 
I F ( 1QUIT.GT.0)  GO  TO  170 
GO  TO  55 
100  CONTINUE 

INCREASE  INTEGRATION  INTERVAL  TO  CALCULATE  FEWER  POINTS. 

QUIT  AFTER  5  TRIES. 

I F ( IFEW.EQ.l )  GO  TO  121 
I F ( I TRY.GT . 5 )  GO  TO  106 
ITRY=I TRY+1 
IMANY=1 

XLIMIT=4.*XLIMIT 
DELW=4./57. 29578 
GO  TO  5 

106  WRITE(6»110)  IDEN 

110  FORMAT ( 16H1 IDENT I FICATION  2A6 , 10X  .  13HPRQGRAM  STOP  / 

160H0THERE  ARE  TOO  MANY  LAMINAR  POINTS  BEING  CALCULATED  -  RN/RB  PRO 
2BABLY  TOO  SMALL  -  ) 

111  I TRBLE=1 
N= J 
N1=N 

CALL  OUTPUT 
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I F ( IMANY.NE.l )  GO  TO  150 
S I GMA=S I GM A/57.295 
I  TRBLE  =  0 
WR I TE ( 6  *  1 13  ) 

113  FORMAT (39HOTHIS  RUN  WILL  BE  REPEATED  WITH  RN/RB=0  ) 

RN  =  0. 

RNRB=0. 

GO  TO  195 

DECREASE  INTEGRATION  INTERVAL  TO  CALCULATE  MORE  POINTS. 

QUIT  AFTER  5  TRIES. 

115  I F  (  I MANY • EQ • 1 )  GO  TO  106 
I F ( I TRY.GT • 5 )  GO  TO  121 
I TRY= I TRY+1 
I  FEW=  1 

XLIMIT=XLIMlT/4. 

DEL W=(4./2.** (FLOAT ( I  TRY ) *2 . ) 1/57.29576 
GO  TO  5 

121  I F ( J.GT.I )  GO  TO  170 
WR I TE ( 6 i 130 )  IDEN 

130  FORMAT ( 16H1IDENTIFICATI0N  2A6  »10X *  1 3HPR0GRAM  SrOP  / 

179H0THERE  ARE  TOO  FEW  LAMINAR  POINTS  BEING  CALCULATED  -  RN/RB  PROB 
2 ABLY  TOO  LARGE  -  ) 

GO  TO  111 

DO  LINEAR  INTERPOLATION  FOR  LAST  POINT. 

170  CF ( J ) =X I  NT ( SMAX  »S ( J— 1 ) *S < J ) »CF ( J-l ) ,CF ( J ) > 

X(  J)=  XINT(SMAX»S( J-l ) , S ( J ) . X ( J- 1 ) * X ( J  )  ) 

EM1( J)=XINT(SMAX,S( J-l) .S( J) ,EM1( J-l) ,EMl( J) ) 

T1 ( J)=XINT(SMAX,S( J-l ) » S ( J ) » T 1 ( J-l )  *T1(J)  ) 

RC ( J ) =X I  NT ( SMAX  »S ( J~1 ) »S ( J ) »RC ( J-l ) »RC ( J ) ) 

RTHC( J)=XINT(SMAX,S( J-l ) » S ( J ) » R I HC { J— 1 ) » R I HC ( J ) ) 

RTH1(J)=XINT(  SMAX  » S ( J-l ) »  S ( J ) ,RTH1< J-l) .RTH1 ( J) ) 

DELTA ( J)=XINT(SMAX,S( J-l) ,S< J) .DELTA (J-l) » DELTA ( J  )  ) 

RES(  J)=XINT(SMAX,S( J-l) ,S( J) ,RES( J-l! »RES ( J ) ) 

FI ( J)=XINT(SMAX,S( J-l ) »S< J) »F1 (J-l >  * F 1 { J ) > 

CIBARlf J)=XINT(SMAX,S( J-l) »S ( J) ,C IBARl ( J-l ) , C I  BARI ( J ) ) 

S ( J ) =SMAX 
150  N= J 

RETURN 

END 
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SIBFTC  TURB.  DD.DECK 
SUBROUTINE  TURB 

C  THIS  ROUTINE  CALCULATES  LOCAL  PROPERTIES  FOR  THE  TURBULENT  PORTION 
COMMON  DELM  ,OELRC  ,OELW ,  1QU 1  T  ,  I  SK I P  , I  I  RbLE  ,  I  TRY  ,K » KODE ,KEY  . 

*  LINES  .MODE  »N»N1,XLIMIT 

COMMON  EMC , EM I N ,GAM , l DEN ( 2 ) , PC  ,  P I N , RB , RLTR  ,RNRB , S I 6MA » T I N . T W 
COMMON  CDo.CDF  ,COFL,COFT ,' UF T TL ,CDN »CDP »CDPL .COPT ,CDPTTL»CDTTL» 
*CDTTLB , UCUF I P  »DCDF TC  » DCDP 1 P ,C I  BARA , TC 
COMMON  CF( 500 ) »CIbARl( 500) , DELTA! 500) »EM1 ( 500  )  »RC ( 500)  ♦ 

*RES( 500) ,RTHC( 500 ) , RTH1 ( 5 00 ) ♦ S ( 500 ) »  T 1 ( 500 ) *X ( 500 ) 

COMMON  CONST ,CPC , C0SWS2 , DELS , D IN , DOVTH * DTHDM1 »GM1 ,GP1» 

*  OMEGA, PCPIN.P  I  ,  RBC,RCTR,RN,S1NWS2,SMAX,SMAXT» 

*ST ,  T  ANW2  »TANwS2»TPTl,TTl,XFl,XF2»  XMAX , XMPM1 , XMU1MW, XM1MC 
COMMON  FI (500) .OMEGAS ( 500) ,XRRR( 750) »YRRR(750) 

EQUIVALENCE  (PC.Pl) , ( PCP I N , P 1 P I N ) 

C 

C  SET  N  USED  IN  EQ.  (12)  AND  (13),  REF.  6  EQUAL  TO  7. 

DATA  EN/7 • / 

FUNCTION  STATEMENT  FOR  LINERAR  INTERPOLATION. 

XINT(XPT,X1 ,X2,Y1 , Y2 ) = Y1+ ( Y2-Y 1 ) * ( XPT-X1 ) / ( X2-X1 ) 

J=N 

DELM=.05 

USING  EQ.  (21),  REF.  6,  CALCULATE  R  THETA(C)  At  A  NEW  POINT  USING 
VALUES  AT  A  PREVIOUS  POINT  FOR  C(F),  F(l),  R  THETA(C),  R(C)  AND 
A  PREDETERMINED  DELTA  R(C). 

10  RTHC( J+1)=RTHC< J)+(CF( J)/(2.*< 1.+F1 (U, \) -RTrtC ( J ) /RC ( J) ) *DELRC 

INCREMENT  R ( C )  BY  DELTA  R(C). 

RC( J+1)=R:< Jl+DELRC 
I  SKI P=0 

START  ITERATION  FOR  LOCAL  MACH  NUMBER  M(l)  AND  OMEGA(S)  USING  EQS. 
(22)  AND  (17),  REF.  6. 

INCREMENT  M(l)  AT  PREVIOUS  POINT  BY  DELTA  M ( 1 )  AND  USE  AS  A  START¬ 
ING  VALUE,  CALL  IT  M ( 1 )  FIRST. 

EM1(J+1)=EM1( Jl+DELM 
K=  1 

30  Tl(J+l)=TIN*(2«+GMl*EMIN**2)/(2.+GMl*EM) { J+l)**2) 

TT1=T1( J+l)*( l.+GMl*EMl(J+l)**2/2. ) 

CALCULATE, OMEGA(S)  FROM  EQ.  (22),  REF.  6. 

DTH=(EN+l.)*(EN+2.)/<EN+2.*(TTl-TW)/TTl)-(EN+2.)*(l.-(TTl-TW) 
1/TT1)/EN+GMi*EM1 ( J+l)**2*( .34+.18M 1 . - ( TT 1-TW > / TT 1 ) )/2. 

I F ( ISKIP.LT. 0)  GO  TO  70 

TEMP  1  =  1.1 2  5*EM I N*SQRT ( ( 2 . +GM1*EMI N**2 ) / ( 2 ,+GMl*EMl ( J+l ) **2 ) )/ 
1(SIN(SIGMA)*P1PIN*EM1(J+1) ) 

TEMP2  =  TEMP1/ ( DTH*RC ( J+ 1 ) *R  THC  <  J+ 1 ) /RbC**2 ) 

TEMP3=SQRT ( TEMP2+TANW2 ) 

OMEGAS ( J+l )=ATAN( TEMP3) 

PLUG  OMEGA ( S )  INTO  EQ.  (17)  REF.  6,  AND  CALCULATE  A  NEW  Mil).  CALL 
IT  M ( 1 i  SECOND 
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SI NWS2=S IN (OMEGAS < J+l) )**2 
C0SWS2=C0S ( OMEGAS ( J+l ) )**2 
TEMPl=GPl/(2. #GAM*EMI N*#2*S I NWS2-GM1 ) 

I F ( TEMPI . 3T.0.  )  GO  TO  40 
EM1( J+1)=2.*EM1(J)-EM1( J- l > 

OMEGAS! J+l ) =2.*0MEGAS ( J ) -OMEGAS ( J-l ) 

ISKIP=-1 
GO  TO  30 

40  TEMP2=TEMP1**< l./GAM) 

TEMP3=GP1*EMIN**2*SINWS2*( 2 .+GM1*EMI N**2 ) / ( 2 . * ( 2 . +GM1*EMI N**2*S INW 
1S2 )  ) 

TEMP 4= ( TEMP2*TEMP3/P1PI N*# ( GM1 /GAM )-l.)*(2./GMl) 

EM1TMP=SQRT ( TEMP4 ) 

DEL1=EM1TMP-EM1( J+l) 

COMPARE  THE  M<1)  FIRST  AND  Mil)  SECOND. 

1)  IF  AbSlM(l)  FIRST-M(i)  SECOND)  LESS  THAN  .001  ACCEPT  THE  VALUE 

2)  IF  GREATER  THAN  .001  SET  Mil)  FIRST  =(M(1)  FIRST  +M(1)  SECONDJ/2 
AND  TRY  AGAIN. 

3)  IF  ITERATION  FAILS  AFTER  100  PASSES.  PRINT  WARNING,  ACCEPT 
CURRENT  VALUE  AND  CONTINUE  TO  NEXT  POINT. 

I F (  ABSIDELl ) .LE.l.E-3)  GO  TO  60 
EMI  I J+l )  =  { EMI { J+l ) +EM1TMP ) /2. 

K=K+1 

IF(K.LE.IOO)  GO  TO  30 
WR I TE { 6 »  50 )  IDEN.X(J) 

50  FORMAT! 16H0 1  DENT  I F I  CAT  I  ON  2A6/91H0wARN ING  -  SUBPROGRAM  TURB  -  ITER 
*AT I  ON  FOR  OMEGAS  AND  Ml  FAILED  AFTER  100  ITERATIONS  AT  X=  1PE11.2) 
60  EMI! J+1)=EM1TMP 
I5KIP=-1 
GO  TO  30 
70  J=J+1 

CALCULATE  T ( PR  I  ME ) / T ( 1 ) 

TPT1= ( 1 .+. 110*EM1 ( J)**2 )*( l.-( .276+.0148*EMl ( J) )*  < TT1-TW) /TT1-.449 
1*< (TT1-TW)/TT1)**3) 

CALCULATE  MU ( 1 ) /MU ! WALL ) 

XMU1MW=(T1( J)/TW)**1.5*(TW/T1 ( J ) +198 . 7/T1 { J) i / ( 1 . +198. 7/T 1! J ) ) 
CALCULATE  R  THETA ( 1 )  EQ.  (22),  REF.  5 

RTH1! J)=RTHC( J ) *tMl ( J ) * ( T1 ( J ) +198 . 6 ) * ( TC/T 1 ( J) ) **2/ ( EMC* ( TC+198 . 6 
1  )  ) 

CALCULATE  THE  DERIVATIVE  D ( DELTA /THETA-H ) /DM1  FROM  EQ.  (14),  REF. 

6.  NOTE  -  EQ.  IS  INCORRECT  IN  REPORT.  (TT1-TW)/TT1  SHOULD  READ 
TW/TT1. 

DTHDM1=GM1*EM1 ( J ) * (  .34+ , 18  +  TW/ TT 1 ) 

C  CALCULATE  F ( 1 )  EQ. ( 19 ) »  REF.  6. 

TEMPI = ( 2 .+GM1*£M1 lJ)**2)-DTH 

TEMP2=1.+GM1*EM1( J) **2+EMl ( J)*!2.+GM1*EM1 ( J)**2)/ !2.*DTH)*DTHDM1 
TEMP3=EM1 ( J)**2*TANWS2/COSWS2*(2.*GAM*cMIN#*2*SINWS2-GM1)*(2.+GM1 
1*EMIN**2*SINWS2)/ ( 2 .* ( TANWS2-TANW2 ) * ( EM  I N*#2*S I NWS2-1  )**2) 

FI ( J ) =-TEMPl / ( TEMP2-TEMP3 ) 
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C 

C  CALCULATE  X  FROM  R(C)  AND  CONE  PROPERTIES. 

X(J)=RC(J)/(222.*PC*EMC*1.E6*(TC+198»6)/TC**2) 

C 

C  CALCULATE  S  USING  X  AND  BODY  GEOMETRY. 

S(J)=RN*(PI/2.-SI GMA ) +X ( J ) -RN / T AN  <  S I GMA ) 

C 

C  CALCULATE  REYNOLDS  NUMBER  BASED  ON  LOCAL  PROPERTIES  AND  S. 
RES!J)=222.*P1*EM1(J)*S!J!*<T1<j)+198w6)*1.E6/T1!J)**2 

CALCULATE  DELTA/THETA  EQ.  (12).  REF.  6. 

DOVTH= (EN+l,)#(EN+2.)/(EN+2.*<TTl-TW)/TTl)+GMl*EMl(J)**2*(1.34 
1+1.38*!1.-(TT1-TW)/TT1) 1/2. 

CALCULATE  DELTA 

DELTA(J)OOVTH*RTHC(J)*X(J)/RC(J) 

CALL  SUBROUTINE  TO  CALCULATE  LOCAL  C { F  » . 

CALL  CFCALC(J) 

DELM=EM1(J)-EM1(J-1) 

I F ( S ( Jl-SMAX)  10,90.80 

INTERPOLATE  FOR  LAST  TURBULENT  POINT  ON  BODY. 

80  CF ( J ) =X I  NT ( SMAX  »S ( J-l ) , S < J ) ,CF < J-l )  ,CF ( J ) ) 

X( J)=  XINT(SMAX,S( J-l ) ,S( J) ,X< J-l) ,X( J) ) 

EM1< J)=XINT(SMAX»S( J-l) »S< J)»EMl(J-l) ,EMl(J) ) 

T1 ( J)=XINT(SMAX,S( J-l ) ,S( J) ,T1 ( J-l ) ,T1( J) ) 

RC ( J) =X I  NT ( SMAX , S (J-l) ,S(J) »RC(J-1)  ,RC(J) ) 

RTHC( J)=XINT(SMAX,S( J-l) ,S( J) , RTHC < J-l ) , RTHC ( J ) ) 

RTH1( J)=XINT(SMAX,S(J-1) ,S< J> ,RTH1<J-1)  ,RTH1( J)  > 

DELTA! J)=XINT<SMAX,S( J-l) ,S(J) .DELTA ( J-l ) .DELTA < J  )  ) 

RES! J)=XINT( SMAX, S( J-l) ,S(J) , RES (J-l)  .RES! J)  ) 

S( J)=SMAX 
90  N  =  J 
RETURN 
END 
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$  IBF.TC  SUMUP.  DD 

SUBROUTINE  SUMUP 
DIMENSION  SUM ( 6 ) 

COMMON  DELM.DELRC.DELW, IQU I T , I SK I P , I TRbLE , I  TRY ,K » KODE »KE Y » 

*  LINES. M0DE»N,N1,XLIMIT 

COMMON  EMC,  EM  IN,  GAM,  1DEN(2  )  , PC  ,P  I N ,  Rt)  ,RLTR  »RNRB ,S  I GMA,  T  IN  »TW 
COMMON  CDb,CDF,CUFL,CDFT,CuFTTL,CuN,CDP,Ci>PL«CDPT,CDPTTL,CD7TL, 
*CDTTLB  »DCDF  I  P  ,uCDF  TC ,  DCUP I  P  ,C  I  tiARA ,  TC 
COMMON  CF< 500), Cl  BARI (500) .DELTA (500) ,EM1(500) »RC(500) » 
*RES(500) ,RTHC( 500) »RTH1 (500) ,S(500) , Tl ( 500 ) *X ( 500 ) 

COMMON  CONST, CPC, C0SWS2 , DELS, 0 1 N , D0VTH.DTHDM1 .GM1.GP1. 

*  OMEGA, PCPIN, PI ,  R8C , RCTR  »RN ,S INWS2 ,SMAX ,SMAXT » 

*ST » T  ANW2 »TANWS2»TPT1»TT1»XF1,XF2» XMAX ,XMPM1 ,XMU1MW,XM1MC 

COMMON  FI (500) .OMEGAS ( 500 ) ,XRRR< 750 ) »YRRR(750) 

EQUIVALENCE  (PC,P1) , ( PCP I N , PI P IN ) 

DO  10  1=1,6 
10  SUM! I ) =0. 

USE  RECTANGULAR  INTEGRATION  FOR  THE  FOLLOWING  INTEGRALS. 

DO  20  I =N1 ,N 

TERM1=RN+(S( I )-ST)*SlN( SIGMA ) 

TERM=P1P IN* ( EMI ( I ) /EMIN )**2*T£RM1*C0S (SIGMA ) *DELS 

EVALUATE  INTEGRAL  FOR  PRESSURE  DRAG  COEFFICIENT. 

SUM( 1)=SUM(1 )+CPC*TERMl*SIN (SIGMA )*DLLS 

EVALUATE  INTEGRAL  FOR  SKIN  FRICTION  DRAG  COEFFICIENT. 

SUM ( 2 ) =SUM (2)+TERM*CF( I ) 

I F ( MODE .EQ.2.AND.N1.NE.1)  GO  TO  20 
I F i MODE . EO. 3 )  GO  TO  20 

FOR  LAMINAR  FLOW  ONLY. 

DE=.968*TW/(EM1( I )**2*T1( I ) )+.145*GMl 

EVALUATE  INTEGRAL  FOR  AVERAGE  LOCAL  CHI  BAR  PARAMETER. 

SUM( 3 ) =SUM (3)+TERMl*CI BAR  1 ( I )*DELS 

EVALUATE  INTEGRAL  FOR  INCREMENTAL  INCREASE  IN  THE  SKIN  FRICTION 
DRAG  COEFFICIENT  DUE  TO  TRANSVERSE  CURVATURE. 

SUM  (  4 )  =SUM  (4)+TERM*(CF(  I  )  *  ( .  5 1 7+ .  91  3*TW/T  1  (  I  )+»121*GMl*EMHI)** 
1*C IBAR1 ( I ) /(EMI ( I )**3*TAN( SIGMA) *  SORT ( 3.  )  )  ) 

EVALUATE  INTEGRAL  FOR  INCREMENTAL  INCREASE  IN  THE  SKIN  FRICTION 
DRAG  COEFFICIENT  DUE  TO  INDUCED  PRESSURE. 

SUM ( 5 ) =SUM (5)+TERM*(CF( I )*(-l,152+.733*TW/Tl( I )+( . 528*GAM-. 126 > 
1*  FMl ( I )**2 ) *OF*XF 1*0 l BAR  1 ( I ) / ( GAM* EMI ( I ) **? ) ) 

EVALUATE  INTEGRAL  FOR  INCREMENTAL  INCREASE  IN  THE  PRESSURE  DRaG 
COEFFICIENT  DUE  TO  INDUCED  PRESSURE. 

SUM(6)=SUM(6)+PCPIN*(XF1*DIN*CIBAR1( I ) +XF2*DI N**2*C IBAR1 ( I ) **2 ) 
1*TERM1*SIN(SIGMA)*DELS 
20  DELS=S ( I+l)-S ' I ) 

IF ( MODE. EQ. 2. AND. Ml .NE. 1 )  GO  TO  30 
I F ( MODE . EQ. 3 )  GO  TO  30 
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CIBARA=2. »SUM(3)/CONST 
DCDFTC=2.*SUM(4)/RB**2 
DCDFIP=2.*SUM< 5)/RB**2 
DCDP1P=4.*SUM(6)/ <GAM*EMIN**2*Rb**2) 
30  CDP=?.*SUM(1)/RB**2 
CDF=2 .*SUM ( 2 ) /RB**2 
RETURN 
END 


$  I BFTC  CFALT.  DD 

SUBROUTINE  CFALT ( J ) 

THIS  ROUTINE  CALCULATES  A  LOCAL  SKIN  FRICTION  COEFFICIENT  FOR 
TURBULENT  FLOW  WHENEVER  THE  I TtRATION  IN  SUuROUTINE  CFCALC  FAILS 
TO  CONVERGE.  THE  METHOD  USED  IS  CONSIDERED  SOMEWHAT  LESS  ACCURATE 
THAN  THE  ITERATIVE  METHOD. 

COMMON  DELM.DELRC.DELW, IQUIT, I SK I P  ,  I TRBLE , I  TRY ,K » KODE , KEY , 

*  LINES, MODE, N,N1,XLIMIT 

COMMON  EMC, EMIN, GAM, I  DENI  2) ,PC »P I N ,Ro  ,RLTR ,RNRB ,S IGMA, T I N » TW 
COMMON  CDB,CDF,CDFL,CDFT,CDFTTL,CDN,CDP,CDPL,CDPT,CDPTTL,CDTTL, 
*CDTTLB,DCDFIP,OCDFTC,DCDPIP,CItiARA,TC 
COMMON  CF ( 500 ) ,CIBAR1(500) ,D£LfA(500>  , EMI ( 500 5 ,RC < 500 ) , 

*RES < 500 ) ,RTHC< 500) ,RTH1 (5005 ,S( 500 > , T 1 ( 500 ) j X ( 500 ) 

COMMON  CONST, CPC, C0SWS2, DELS, D IN, DGVTH,DTHDM1,GM1,GP1, 

*  OMEGA, PCPIN, PI ,  RBC ,RCTR  »RN  ,S INWS2 ,SMAX ,SMAXT , 

*ST  ,TANW2 , TANWS2 »TPT1,TT1»XF1,XF2» XMAX , XMPM1 »XMU1MW»  XM1MC 

COMMON  FI ( 500 ) .OMEGAS ( 500 ) ,XRRR ( 750 ) »YRRR(750) 

EQUIVALENCE  ( PC ,P1 ) , ( PCP I N ,P1 P IN ) 

TERM1=ALOG10(2.*XMU1MW*RTH1 ( J) ) 

TERM2=( .8686+TERMl)*TERMl*TPTl 
CF(J)=. 05856/TERM 2 


I F ( KEY «  EQ. 10 )  WRI TE ( 6 , 5 ) 

5  FORMAT  (1HU 
KEY=20 
WR I TE ( 6 , 10 ) 

10  FORMAT (40H  ALTERNATE  METHOD  USED  TO  CALCULATE  C(F)  ) 
RETURN 
END 
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$  I BFTC  CF.  DD 

SUBROUTINE  CFCALC(J) 

C  THIS  ROUTINE  CALCULATES  LOCAL  SKIN  FRICTION  COEFFICIENT  FOR 
C  TURBULENT  FLOW. 

COMMON  DELM.DELRC.DELW, IQU I T , I SK 1 P , I TRBLE , I  TRY ,K » KODE , KEY  * 

*  LINES, MODE, N,N1,XLIMIT 

COMMON  EMC,  EMIN, GAM,  I  DEN  (  2)  ,PC , P I N , Rb ,RLTR »RNRb ,S I GMA, T I N , TW 
COMMON  CDB,CDF,CDFL»CDFT,GuFTTL,CuN,LDP,CDPL,CDPT ,CUPTTL»CDTTL, 
#CDTTLb»DCDFIP,DCDFTC»DCDPIP*CIBARA»TC 
COMMON  CF ( 500 ) ,C I  BARI (500) .DELTA (  500) ,EM1 ( 500 ) , RC ( 500 > » 

*RES ( 500 ) ,RTHC( 500 ) ,RTH1 (500) ,S( 500)  , T 1 ( 500 ) , X ( 500  ) 

COMMON  CONST, CPC, C0SWS2, DELS, D I N,D0VTH,DTHDM1 ,GM1  ,GP1, 

*  OMEGA ,PC°IN  ,P I ,  RBC,RCTR,RN,SINWS2,SMAX,SMAXT, 

*ST , TANW2 , TANWS2,TPT1 ,TT1,XF1,XF2,XMAX,XMPM1,XMU1MW,XM1MC 

COMMON  FI ( 500) , OMEGAS ( 500) ,XRRR( 750) ,YRRR(750) 

EQUIVALENCE  (PC, PI ) , ( PCP I N , P 1 P IN ) 

SET  PRANDTL  NUMBER  USED  IN  EQ.  (6),  REF.  6  EQUAL  TO  .75 
DATA  PR/0.75/ 

START  ITERATION  FOR  C(BIG  F)  AND  ETAtBIG  F) 

IF  THIS  IS  THE  FIRST  PASS  THRU,  SET  tTA ( D IG  F)  =  .5 
AS  A  FIRST  APPROXIMATION.  CALL  IT  ETAF1 
I F ( KEY . EQ.O )  ETAF1= . 50 
KEY=10 

CALCULATE  BETA  FROM  EQ.  (2),  REF.  6. 

BETA=GM1*EM1( J ) **2/ ( 2 .+GM1*EM1 (J)**2> 

CALCULATE  RECOVERY  TEMPERATURE  TE1  EQ.  (6),  REF.  6 
TE 1 = ( l.+PR**< 1./3. )*GM1*EM1( J)**2/2. )*T1( J) 

CALCULATE  SUBLAYER  PARAMtTER  S  EQ.  (5),  REF.  6 
SS =11.5+6. 6* (TE1-TW)/TE1 

CALCULATE  PHI  FROM  DEFINITION  PAGE  4,  REF.  6 
PH  I  =SQRT  (4.*bETA*(l.-(TTl-Tvii)/TTl)  +  !  (  TT1-TW )  /  T  T 1 )  **2  ) 

CALCULATE  ETA  FROM  DEFI NATION  PAGE  3,  REF.  6 
ETA=(2.*BETA-(TT1-TW)/TT1) /PHI 
KOUNT  =  1 

10  IF(ABS(ETAF1) .LT.l. )  GO  TO  15 

IF  ITERATION  HAS  FAILED,  USE  ALTERNATE  METHOD  TO  CALCULATE  LOCAL 
SKIN  FRICTION  COEFFICIENT. 

11  CALL  CF  ALT ( J ) 

ET  AF1= . 5 
GO  TO  70 

CALCULATE  C(BIG  F)  EQ.  10,  REF.  6 

15  CBF= ( ( ,242*(ARSIN(ETA) -ARS I N ( ETAF 1 ) ) > / ( -1 . 968+ALOG10 ( 23 .*XMU1MW* 
1RTHK J)/S  3) ) )**2/(BETA*TTl/Tl  (J) ) 

PLUG  CALCULATED  VALUE  FOR  C(bIG  F)  INTO  DEFINITION  FOR  ETAtBIG  F) 
PAGE  3,  REF.  6.  CALL  IT  ETAF2 
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ETAF2=(2.*bETA*SS*SQRT(TW*CbF/(Tl{ J)*2.) ) - ( TT 1-TW ) / TT1 ) /PHI 

COMPARE  ETAF2  WITH  ETAF1 

IS  ABSOLUTE  DIFFERENCE  LESS  THAN  .00001 

IF  YES  STOP  ITERATION 

IF  NO  SET  ETAFl=(ETAF2+ETAFl)/2.  AND  CONTINUE 
IF(ABS(ETAF2-ETAFl).LE.l.E-5)  GO  TO  50 
ETAFl=(ETAF2+ETAFl)/2. 

IF ( KOUNT »GT . 1000 )  GO  TO  20 

KOUNT=KOUNT+l 

GO  TO  10 

20  WR I TE ( 6 * 30 )  IDEN.X(J) 

30  FORMAT! 16H0 1  DENT  I F I  CAT  I  ON  2A6/91H  WARNING  -  SUBPROGRAM  CFCALC  -  IT 
DERATION  FOR  ETAF  AND  CF  FAILED  AFTER  1000  ITERATION  AT  X=  1PE11.2) 
50  XRAY=1.-ETAF2**2 

IF (XRAY.GE.O. )  GO  TO  60 
GO  TO  11 

USING  CALCULATED  VALUES  FOR  ETAlblG  F)  AND  CtbIG  F)  SOLVE  EQ.  (91* 
REF.  6  FOR  LOCAL  SKIN  FRICTION  C ( SMALL  F). 

NOTE  -  EQ.  (9)  REF.  6  IS  INCORRECTLY  WRITTEN  IN  THE  REPORT,  THE 
THIRD  TERM  IN  THE  NUMERATOR  AND  DENOMINATOR  SHOULD  BE  MULTIPLIED 
BY  SUBLAYER  PARAMETER  S. 

60  TERM1=-1.968+ALOG10(23.*XMU1MW*RTH1 ( J ) /SS )  +  .242*SS*SQRT ( 2 .*BETA* 
1TW/TT1)/(?HI*SQRT(XRAY) ) 

TERM2=TERM1+1. 968-1. 099 
CF(J)=CBF*(TERM1/TERM2) 

ETAF1=ETAF2 
70  RETURN 
END 
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SUBROUTINE  OUTPUT 

COMMON  DELM  ,DELRC ,DELW, lUU I T  » I SA I P » I TRuL t , I  TRY ,K , KOUt ,KEY  . 

*  LINES, MOuE,N,Nl,XLIMIT 

COMMON  EMC,  EMIN, CAM,  IDEM  2  )  ,PC  ,P  I , Ko ,RLTK  ,Ki'(Ro , S  IGi-.A, T  1 N , T W 
COMMON  CDtS.COF  ,  CDh  L  »CUFT  ,  Cuf- T  TL  ,  CUn  ,  CUP  ,CU>PL  ,CuP1  ,CwPTTL»CUTTL, 
#CDTTLB,DCDFIP,OCU)-TC,DCUPIP,C  IdmRa.TC 
COMMON  CF( 500) ,CI BARI { 500 ) .DELTA ( 500) ,CMI ( 5C0) ,RC(500) , 

*RES (  500 ) ,RTHC( 500 ) ,RTH1 (500) ,S(500! , T 1 ( 500 ) , X ( 500  ) 

COMMON  CONST »CPC,C0SWS2, DELS, DIN,U0VTH,UTI1DM1  ,GM1  ,GP1, 

*  OMEGA, PCPIN, PI  ,  RBC,KCTR,RiN,SIi'U\'S2,5MAX,SMAXT, 

*ST ,  T  ANW2  » TANWS2  ,TPT1,TT1,XF1,XF2,  Xi'iAX  ,XMPM1 ,  Xi'.UlMW,  XhlMC 
COMMON  FI  (  500)  .OMEGAS  (  500)  ,XRRR(  750)  ,YRRR(  7 1» 0  ) 

EQUIVALENCE  (PC, PI) , ( PCPI N ,P1P IN) 

AX  I  AL  =  RN*  ( 1  •— S  IN  (  S I  G*’iA )  )  +  (  Ki-Rn*S  I N  ( P  I  / 2 . -S I Gi'.A )  )  *CGS (  S I GMA  ) 

1  /  S I N  (  S I GMA ) 

RLIN=222.*PIN*cMIN*AXIAL*l.t6*(T I N+ 1 98 . 6 ) / T I N**2 
SIGMA=SIGMA*57. 29578 
L I NES=0 

GO  TO  ( 10,30,40) , MODE 
10  WR I TE ( 6 , 20 )  IDEN 

20  FORMAT(16H1IOEimTIFICATIO.')  2A6,12X,1hh  LAMINaR  UNLY  ) 

GO  TO  50 

30  WR I TE ( 6 , 35 )  IDEN 

35  FORMAT! 1 6H1 1  DENT  I F I  CAT  I  ON  2A6,12X,19H  LAMINAR-TURBULENT  ) 

GO  TO  50 

40  WR I TE ( 6 , 45 )  IDEN 

45  F0R1<1AT(16H1IDEi'iTIFICaTI0N  2A6,12X,11H  TURBULENT  ) 

50  WRITE  (6,60)  P I N ,  T  IN  ,  tM  I  N,  PC ,  T  C ,  L.-.C ,  S  IG.-.a  ,  T  *  ,  RiYRtJ ,  Ro  ,RL  IN  ,RLTR 
60  FORMAT ( 1H0  12X  9HP  ( I NF ) =  1PE 13 . 6 , 7X , 9HT  ( IMF ) =  1PL13.6,7X»9HM  (IN 
IF ) =  1PE13.6/1H0  1 7X  4HPC=  1PE13.6,12X,4HTC=  1PL1 3 .6 , 12X, 4HMC  =  1PE1 
23.6/1H0,  14X,7HSIGMA=  1PE13.6 ,6X  ,  lOHT  (v;ALL)=  1PE 13 .6 , 9X , 7HRN/Rb= 
31PE13.6/1H0,17X,4HRB=  1PE13.6 ,6X , lOHRL  { I NF ) =  1PE13.6, 

4  7X.9HRL  ( TR ) =  1PE13.6) 

L I NES=L I NES+9 
I F ( ITRBLE.GT.O)  60  TO  119 
GO  TO  (65,65,100)  ,  MODE 

65  WR I TE ( 6 , 70 )  DCDP I P , DCDFTC , uCDF I P ,C I ohRa ,CoKL ,CDPL 

70  FORMAT (1H0»7X,14HD  CDP(1P)(L)=  lPtl3.6,2X, 14HU  CDF(TC)(L)=  1Pc13.6 
1 ,2X , 14H  D  CDF ( IP ) ( L ) =  1PE1 3 .6/ 1H0 ,6X , 1 5HCH I  BAR  AVG(L)=  1PE13.6.8X, 
28HCDF ( L ) =  1PE13.6,8X,8HCDP(L)=  1PE13.6) 

L I NES=L I NES+4 
GO  TO  ( 100,80) .MODE 
80  WR I TE ( 6 , 90 )  CDFT , CDPT 
L I NES=L I NES+2 

?0  FORMAT ( 1H042X.3HCDF ( T ) =  1  PL  13 . 6 , 8X , 8HCDP ( T ) =  1PE13.6) 

100  WRI  TE(6,ll0)CDN,CDrTTL,CuPTTL,CUB,CuTTL,CuTTLo 

110  FORMAT ( 1  HO , 1  OX  » 1 1HCD  (NOSE)=  1PE 1 3 . 6 , 3X , 1 3HCDF  ( TOTAL ) =  1PE13.6, 
13X , 13HCDP  ( TOTAL) =  1PL13.6/1H0,2X,19HCU  (  BASE )( VACUUM) =  lPbl3.6»4x 

2  12HCD  ( TOTAL ) =  1PE13 .6  »4X , 12HC0  ( TOTAL ) =  1PE13.6/1H  41X .48HEXCLUD 

3ING  BASE  DRAG  INCLUDING  BASE  DPAG/lHO) 

L I  NES  =  L  I  NF.S+7 

GO  TO  ( 119,119,111 ) ,  MODE 

111  WR I TE ( 6 , 1 1 5 ) 

115  FORMAT ( 1H0  5X,  1HX,  10X  1HS  9X  2HM1  9X  2HCF  11X  2HT1  8X  OHDELTA  8X 
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13HRXC  10X  4HRTHC  9X  4HRTH1  8X  3HRS1  ) 

GO  TO  125 

119  WRITEI6.120) 

120  FORMAT ( 1H0  5X,  1HX  10X  1HS  9 X  2HM1  9X  2HCF  11X  2HT1  8X  5HDELTA  8X 
13HRXC  10X  4HRTHC  9 X  4HRTH1  CX  3HRS1  6X  7HCH I bARl  ) 

125  L I NES=L I NES+2 

IF(NI.EO.I)  Nl=-1 
DO  220  1=1, N 
I F ( I .GT.N1 )  GO  TO  140 

WRITE (6,1 30 )  X(  1  ) »S<  1)  ,EM1 ( I ) ,CF( 1 ) ,T1( I ) , DELTA ( I ) ,RC( I ) »RTHC( I  ) , 
1RTH1 ( I ) »RES ( 1 ) »C I  DARI ( I ) 

130  FORMAT (1P4E11.2,1P7E12»3) 

GO  TO  170 

140  IF( I.GT.(N1+1) )  GO  TO  160 
WRITE(6,150) 

150  FORMAT  (  lOHOTURtiULENT) 

160  WRITE! 6, 130 )X< I ) ,S( I ) ,EM1( 1 ) ,CF( I) s  T1 ( I ) ,DELT«( I ) ,RC( I ) ,RTHC( I ) , 
1RTH1 t I ) ,RES( I ) 

170  L I NES=L I NES+1 

IF(LINES.LT.52)  GO  TO  220 
WRITE(6»180)  IDEN 

180  FORMAT ( 16H 1 1  DENT  I F 1  CAT  I  ON  2A6 ,4X,9HC0NT 1NUED  ) 

GO  TO  ( 190,190,200) , MODE 

190  IFt  I.GT.INl  +  l) )  GO  TO  200 
WR I TE ( 6 , 120 ) 

GO  TO  210 

200  WRITE (6, 115) 

210  UNES=2 

220  CONTINUE 
RETURN 
END 

SDATA 
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